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Abstract 

As predicted by quantum chromodynamics(QCD), around T ~ 190 MeV in the early universe, 
the QCD transition occurs during which the quarks are combined into the massive hadrons. 
This process reduces the effective relativistic degree of freedom, and causes a change in 
the expansion behavior of the universe. Similarly, the e + e~ annihilation occurred around 
T ~ 0.5 Mev has the same kind of effect. Besides, the dark energy also drives the present 
stage accelerating expansion. We study these combined effects on the relic gravitational waves 
(RGWs). In our treatment, the QCD transition and the e + e~ annihilation, each is respectively 
represented by a short period of expansion inserted into the radiation era. Incorporating these 
effects, the equation of RGWs is analytically solved for a spatially flat universe, evolving from 
the inflation up to the current acceleration, and the spectrum of RGWs is obtained, covering 
the whole range of frequency > 10~ 19 Hz. It is found that the QCD transition causes a 
reduction of the amplitude of RGWs by ~ 20% in the range > 10~ 9 Hz, and the e + e~ 
annihilation causes a reduction ~ 10% in the range > 10~ 12 Hz. In the presence of the dark 
energy, the combination of the QCD transition and the e + e~ annihilation, causes a larger 
reduction of the amplitude by ~ 30% for the range > 10~ 9 Hz, which covers the bands of 
operation of LIGO and LISA. By analysis, it is shown that RGWs will be difficult to detect by 
the present LIGO, but can be tested by LISA for certain inflationary models. 

PACS numbers: 04.30.-w, 98.80.-k, 04.62,+v 



1. Introduction 



The existence of a stochastic background of relic gravitational waves (RGWs) is generally predicted in 
inflationary models [1, 2, 3]. Since the relic gravitons decoupled much earlier than the cosmic microwave 
background (CMB) photons, the detections of RGWs would open a new window to the very early universe. 
Unlike gravitational radiations from finite stellar objects, RGWs exist everywhere and anytime, and have a 
wide spreading spectrum, serving as one of the major scientific goals of the laser interferometers gravity- 
wave detections, including the current LIGO [4], VIRGO [5], GEO600 [6], TAMA [7], AIGO [8], and 
the future LISA [9], ASTROD [10], BBO [11], and DECIGO [12]. Moreover, along with the density 
perturbations, RGWs also contribute to the CMB anisotropies and polarizations [13, 14, 15, 16, 17]. For 
instance, the B-polarization of CMB on very large scales can only be generated by RGWs. Through the 
detections of CMB polarizations one may have a chance to obtain the direct evidence of GWs for the first 
time [18, 19, 20, 21, 22]. Therefore, the precise information of RGWs is much desired. 

The spectrum of RGWs depends on the following factors. First, it depends sensitively on the specific 
inflationary models [1, 3, 23]. After being generated, it will be altered by the subsequent stages of 
expansion of the universe, among which notably is the currently accelerating expansion of the universe 
[3, 24]. Finally, it would be further modified by other physical processes occurred in the early universe. 
For example, the neutrino free-streaming [25] has been shown to affect the spectrum of RGWs in the 
very low frequency range (10~ 16 ~ 10~ 10 ) Hz [26]. Although this would modify the contribution to 
CMB polarizations, but would not affect the detections of LIGO [4] and LISA [9] that operate effectively 
at higher frequencies, (10 1 ~ 10 3 ) Hz and (10~ 7 ~ 10°) Hz, respectively. Other important physical 
processes include the QCD transition around a temperature T ~ 190 MeV [27, 28], and the e + e _ 
annihilation around T ~ 0.5 MeV. These two processes will be studied in this work. Leaving aside the 
detail of the QCD transition, which is notoriously complicated and still under study, we are concerned with 
only the thermodynamic property that the transition brings to the system, i.e., a significant drop in the 
relativistic degrees of freedom. Consequently, the constituent components in the stress tensor T^ u also 
change in the Friedmann's equation. Thereby the expansion of the Universe in terms of the scale factor 
a(r) will subsequently change, and the spectrum of RGWs would be modified. The e + e~ annihilation 
has the similar effect. Schwarz [29] studied the QCD transition and the e + e~ annihilation and estimated 
the reduction of the energy density spectrum of RGWs. Watanabe and Komatsu [30] investigated these 
effects with more details, and gave a numerical solution of the energy density spectrum of RGWs. In these 
works, the important effect of the accelerating expansion of the present universe [3] has not been included, 
neither the effects of inflation and reheating. Moreover, to show the prospects of detecting RGWs at the 
laser interferometers gravity-wave detection, an explicit demonstration of the current spectrum of RGWs 
itself and its direct comparison with the sensitivity curve of gravity-wave detection are also needed. 

In this paper, by extending our previous analytical calculation of the spectrum of RGWs, we explore the 
consequences caused by the QCD transition and the e + e~ annihilation, and, at the same time, take into 
account of the accelerating expansion driven by the dark energy [3]. Beside the energy density spectrum 
£l g (y), we demonstrate explicitly the spectrum h(u,Tu) of RGWs itself, and compare it directly with the 
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a comprehensive compilation, by using the set of parameters (3, (3 S , n, v, 7 and r, respectively, such 
important cosmological elements have been explicitly parameterized, as the inflation, the reheating, the 
QCD transition, the e + e~ annihilation, the dark energy and the tensor/scalar ratio. This will considerably 
facilitate further studies on the RGWs and the relevant physical processes. For instance, it can be easily 
used in the calculation of CMB anisotropies and polarizations generated by RGWs [16, 21]. 

The organization of this paper is as follows. In section 2, from the inflation up to the acceleration, 
the scale factor a(r) is specified by the continuity conditions for the subsequential stages of expansion. 
The periods of the QCD transition and the e + e~ annihilation are also modelled as having a scale factor 
a(r) of power-law form, which is inserted into the radiation era. In section 3, the analytical solution of 
RGWs in terms of Bessel's functions is determined with the coefficients being fixed by continuity condition 
joining two consecutive expansion stages. The impact of the dark energy on the expansion is emphasized. 
In section 4, we present the resulting spectra of RGWs and discuss the modifications caused by the QCD 
transition, the e + e~ annihilation, and other effects. Appendix A supplies the details of our treatment of 
the period of QCD transition and e + e~ annihilation by modelling the corresponding scale factor. Appendix 
B gives an interpretation of the modifications on RGWs due to the QCD transition. In this paper the unit 
with c = h = k,B = 1 is used. 



2. Expansion history of the universe 

From the inflationary up to the currently accelerating stage, the expansion of the universe can be 
described by the spatially flat (17a + Q m + J7 r = 1) Robertson-Walker spacetime with a metric 

ds 2 = a 2 (T)[-dT 2 + Sijdx^xj], (1) 

where r is the conformal time. The scale factor a(r) for the successive stages can be approximately 
described by the following forms [31]: 
The inflationary stage: 

o(r) = Iq\t\ 1+(S , -00 <t<ti, (2) 

where 1 + (3 < 0, and t\ < 0. This generic form of scale factor is a simple modelling of inflationary 
expansion, and the index j3 is a parameter. The special case of j3 = — 2 is the de Sitter expansion of 
inflation. If the inflationary expansion is driven by a scalar field, then the index (3 is related to the so-called 
slow-roll parameters, rj and e [32], as (3 = — 2 + (rj — 3e) . In this class of scalar inflationary models one 
usually has (3 < —2. Besides, the observational results of WMAP also indicate that (3 should be slightly 
smaller than —2 [18]. But, for demonstration purpose, we allow the parameter (3 to take values > —2 to 
examine the possibility of detection by GWS detectors. 
The reheating stage: 

o(r) = a z \r - r p | 1+/3 % n < r < r s , (3) 

where t s is the beginning of radiation era and t p < t\. We will mostly take the model parameter 
(3 3 = —0.3, though other values are also taken to demonstrate the effect of various reheating models. 



The radiation-dominant stage: 

The QCD transition occurs around T ~ 190 Mev for a period, and the process of e + e~ annihilation 
into photons starts around T ~ 0.5 MeV and ends up around T ~ 0.1 MeV [30]. These two periods 
should be included into the radiation stage. Before the QCD transition, one has 

a(r) = a e (r - r e ), t s < r < r q ; (4) 

where r q is the beginning of the QCD transition and r e < t s . During the QCD transition, a(r) is modelled 
by 

a(r) = a„(r - r n ) 1+n , T q < r < t x , (5) 

where t x is the ending of the QCD transition and r n < r q . The power index n in Eq.(5) is a model 
parameter describing the QCD transition. We have found that the spectrum of RGWs does not vary 
considerably for the values of n in the interval (0.7 ~ 4.6). For concreteness we take n = 1.634 in 
calculations. (For more details, see Appendix A). The expansion rate a'(r) around T ~ 190 Mev is 
plotted in Fig.l, showing a jump-up caused by the QCD transition. After the QCD transition and before 
the e + e~ annihilation, one has 

a(r) = a/(r - r/), t x < r < Ty , (6) 

where r y is the beginning of the e + e~ annihilation and t/ < t x . The two slopes a e in Eq.(4) and a/ in 
Eq.(6) are related as a/ ~ 1.2a e by considerations of the details of the QCD transition (See Appendix A). 
Similarly, the period of e + e~ annihilation is modelled by 

a(r) = a„(r - t v ) 1+v , r y < r < r z , (7) 

where r 2 is the ending of the e + e~ annihilation and t v < r y . The power index v in Eq.(7), as a model 
parameter, can be taken as v = 0.063 (see Appendix A). After the e + e~ annihilation, one has 

a(r) = a g (T -T g ), t z <t<t 2 , (8) 

where r g < t z , and T2 is the beginning of the matter era, which can be taken at a redshift z ~ 3454 [18]. 
The two slopes aj in Eq.(6) and a g in Eq.(8) satisfy the relation a g ~ 1.1a/ (See Appendix A). 
The matter-dominant stage: 

a(r) = a m (r - T m ) 2 , t 2 <t<t e , (9) 

where r m < t 2 and te is the beginning of the acceleration era. 
The accelerating stage: 

a(r) = I h \t - r a | -7 , t e <t<t h , (10) 

where Tf/ is the present time and th < r a . The index 7 in Eq.(10) depends on the dark energy J1a- By 
fitting with the numerical solution of the Friedmann equation [3, 26], 

/a'\2 8vrG. 
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where a' = da/dr, one can take 7 ~ 1.05 for Q\ = 0.7, and 7 ~ 1.044 for Da = 0.75. The redshift of 
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constant model with ft^ = 0.72 and VL m = 0.28, it starts at z ~ 0.37, and, in the quantum effective 
Yang-Mills condensate dark energy model, it starts at z ~ 0.5 [33]. 

In the above specifications of a(r), there are nine instances of time, 77, t s , T q , t x , r y , t z , t 2 , te, and th, 
which separate the different stages. Eight of them are determined by how much a(r) increases over each 
stage based on the cosmological considerations. We take the following specifications: Ci = = 300 for 
the reheating stage, ( s = ^| = 2.895 x 10 16 for the first radiation stage, ( q = ^) = 1.34 for the second 
radiation stage, ( x = = 283.582 for the third radiation stage, ( y = = 5 for the fourth radiation 
stage, ( z = ^ = 1.818 x 10 4 for the fifth radiation stage (see Appendix A), C2 = = 3454Q 1 for 
the matter stage, and Qe = °(r^) = (r^) 1 ^ 3 ^ or tne P resent accelerating stage [3, 26]. The remaining 
time instance is fixed by an overall normalization 

\TH ~ Ta\ = 1. (12) 

There are also 22 constants in the expressions of a(r), among which 0, f3 s , n, v and 7 are imposed as the 
model parameters describing the inflation, the reheating, the QCD transition, the e + e~ annihilation and 
the acceleration, respectively. Based on the definition of the expansion rate Hq = ^i\ TH of the present 
universe , one has In = j/Hq. Making use of the continuity conditions of a(r) and of a(r)' at the eight 
given joining points n, t s , T q , t x , t v , t z , T2, te, and th, all parameters are fixed as the following: 

Ta-TE = Q, 
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TE ~T m = - Q, 

7 
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and 
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In the expanding universe, the physical wavelength is related to the comoving wavenumber k by 

A=^, (15) 
and the wavenumber ku corresponding to the present Hubble radius is 

l/-no 

There is another wavenumber involved 

_ 2ira(T E ) k H 

kE = 1 ITT = TX^' (17) 
L/Hq I + Ze 

whose corresponding wavelength at the time te is the Hubble radius I/Hq. In the present universe the 
physical frequency corresponding to a wavenumber k is given by 

v = — = — - = k. (18) 



3. Analytical solution of RGWs 
In the presence of the gravitational waves, the perturbed metric is 

ds 2 = a 2 (T)[-dT 2 + (Sij + hijUiydr'. (19) 
where the tensorial perturbation hij is a 3 x 3 matrix and is taken to be transverse and traceless 

h\ = 0, htjj = 0. (20) 

The wave equation of RGWs is 

f)..(,r^7,f) v h,A = 0. m\ 



We decompose into the Fourier modes of the comoving wave number k and into the polarization state 

Mr.x) = E/ T^^^We**, ( 22 ) 



<7 as 



(27)3 

where /i^*(r) = h k a \r) ensuring that /i^ be real, e?- is the polarization tensor, and a denotes the 

» 



polarization states x,+. Here hij is treated as a classical field. In terms of the mode h k , Eq.(21) 
reduces to 

4 CT) 'V) + 2^4 CT) '(-) + k 2 ht\r) = 0. (23) 

By assumption, for each polarization, x, +, the wave equation is the same and has the same statistical 
properties, so that the super index (a) can be dropped from from now on. Since for all the stages of 
expansion the scale factor is of a power-law form 

o(r) oc r a (24) 

the solution to Eq.(23) is a linear combination of Bessel function J v and Neumann function N v 

h k (T)=A- a \a 1 J a _i_{kT) + a 2 N a _ L {kT)\, (25) 



where the constants a\ and a 2 are completely determined by the continuity of h k and of h' k at the joining 
points ti, t s , T q , t x , T y , t z , t 2 , and TE- 

In particular, we write down explicitly the solution for the inflationary stage since it give the initial 
condition for the spectrum of RGWs, 



h k (r) = A 1 1 \t\ 2 (S [A 1 Ji + p{x) + A 2 J_ ( i + p ) {x)~\, -oo < r < n (26) 
where x = kr and 



ii = — Vv? 6 ^ 72 ' A * = -^r-S e ~ lnm > ( 2? ) 

COS (jTT V 2 COS /57T V 2 

are taken [34], so that the so-called adiabatic vacuum is achieved: lim^oo /ifc(r) oc e _lfcr in the high 
frequency limit [35]. Moreover, the constant Aq in Eq.(26) is independent of A;, whose value is determined 
by the initial amplitude of the spectrum. For kr <S 1 the /c-dependence of %(t) is given by 

h k (r) oc Ji +p (x) oc k^ +l3 . (28) 

As will be seen later, this choice will lead to the required scale-invariant initial spectrum in Eq.(32). 

It should be mentioned that around the temperature T ~ 2 Mev neutrinos decoupled from electrons 
and photons and started free-streaming in space. This will give rise to an anisotropic part of the energy- 
momentum tensor as a source of the equation of RGWs. The previous works have shown that the neutrino 
free-streaming in combination with the dark energy would cause a reduction of the amplitude of RGWs by 
~ 20% in the low frequency range (1CT 16 ~ 1(T 10 ) Hz [25, 26, 30]. Although this modified RGWs will 
contribute to the CMB anisotropies and polarizations on very large scales, the frequency range is outside 
the frequency bands of LIGO and LISA. On the other hand, as will be seen, the impact on RGWs by the 
QCD transition is in the high frequency range v > 10~ 9 Hz. Since these two frequency ranges are not 
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In our previous study [26], we have examined the issue of how the RGWs would be, if our current 
universe were matter-dominated. The amplitudes of RGWs in the ACDM (accelerating) and in the CDM 
universe were compared, and the ratio was found to be hk(T}{)hCDM /hk(Tn)cDM ~ 1-3. Here a similar 
examination is extended to the case including the QCD transition and the e + e~ annihilation. To be specific, 
we assume that both universes have the same initial afa) and a'fo) at the time T2 with z ~ 3454, when 
PA *C Pm = Pr- By numerically solving the Friedmann equation in both models, we plot the scale factor 
a(t) in Fig. 2, showing that the ratio of scale factors at present is cl{th)acdm / 'o,{th)cdm — 1-8. As is 
known [1, 3], for wavelengths shorter than the horizon the modes decay as hk(r) oc l/a(r), so the CDM 
model would predict an amplitude of RGWs higher than the ACDM model. This is indeed confirmed by 
our calculation including the QCD transition and the e + e~ annihilation, and the ratio is 

hk {t~h)cdm / hk {t~h) acdm ~ 1-8. (29) 

Moreover, there are some subtleties with the matter-dominant model, regarding to interpreting the current 
observations. As it stands, the actual universe is ACDM, so the observed Hubble constant is properly 
interpreted as the current expansion rate in the accelerating model, Ho = (a' /a 2 ) T[I . The virtual matter- 
dominant universe would have a smaller rate (a'/a 2 ) TH ~ 0.4 Hq. If the observed Hubble constant Hq 
were regarded as the expansion rate of the virtual matter-dominant universe [30], one would come up with 
an amplitude of hk(Tn) lower by an extra factor ~ 1.8 than it should be. 



4. Spectrum of relic gravitational waves 
The spectrum of RGWs h(k,r) at a time r is defined by the following equation [31]: 

roc fjU 

j o h\k,T)^^(0\h^(x,r)h l3 (x,T)\0}, (30) 
where the right-hand side is the expectation value of the h^hij. Calculation yields the spectrum at present 

h(k,r H ) = -k 3 / 2 \h k (r H )\, (31) 

7T 

where the factor 2 counts for the two independent polarizations. 

One of the most important properties of the inflation is that the initial spectrum of RGWs at the time 
Tj of the horizon-crossing during the inflation is nearly scale-invariant [31]: 

h{k, n ) = A(^-) 2+ \ (32) 

where 2 + /3 ~ 0, and A is a constant to be fixed by the observed CMB anisotropies in practice. The First 
Year WMAP gives the scalar spectral index n s = 0.99±0.04, the Three Year WMAP gives n s = 0.951+g ^g 
[18], while in combination with constraints from SDSS, SNIa, and the galaxy clustering, it would give 
n s = 0.965 ±0.012 (68% CL) [20]. The five-year WMAP data give n s = 0.963lffi [36], and the 
WMAP data combined with Baryon Acoustic Oscillations and Type la supernovae give n s = 0.9601q;oi3 
(95% CL) [37].. From the relation n s = 2(3 + 5 [1, 3], the inflation index (3 = -2.02 for n s = 0.951. As 
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Note that the constant A is directly proportional to i n Eq.(26). Since the observed CMB anisotropies 
[18] is AT/T ~ 0.37 x 1CT 5 at I ~ 10, which corresponds to anisotropies on scales of the Hubble radius 
1/Hq, so, as in Refs.[26], we take the normalization of the spectrum 

h(k E ,T H ) = 0.37 x 10~ 5 r5, (33) 

where kg is defined in Eq.(17), its corresponding physical frequency being v% = kEj^aijjj) = H /(l + 
ze) ~ 1 x 10~ 18 Hz. The tensor/scalar ratio r can be related to the slow roll parameter e in the 
scalar inflationary model as r = 16e [32]. However, the value of r is model-dependent, and frequency- 
dependent [16, 17]. This has long been known to be a notoriously thorny issue [38]. In our treatment, for 
simplicity, r is only taken as a constant parameter for normalization of RGWs. Currently, only observational 
constraints on r have been given. The Three Year WMAP constraint is r < 2.2 (95% CL) evaluated at 
k = 0.002 Mpc" 1 , and the full WMAP constraint is r < 0.55 (95% CL) [21]. The combination from 
such observations, as of SDSS, 3-year WMAP, supernovae SN, and galaxy clustering, gives an upper limit 
r < 0.22 (95% CL) [20]. The five-year WMAP gives a limit r < 0.43 (95%) for power-law models [36], 
and the WMAP data combined with Baryon Acoustic Oscillations and Type la supernovae give r < 0.20 
(95% CL) [37]. For concreteness, we take r = 0.22. 

The spectral energy density £l g (k) of the RGWs is given by 

n g (k) = ^-h\k,T H )(^)\ (34) 

which follows from the definition of the total energy density of RGWs [31] 

^GW=^ = [ hUPPer Slg(k)^, (35) 

Pc J how fc 

where p g = 32^3^.7,0^0 is the energy density of RGWs, and p c = 3Hq/8itG is the critical energy 
density. The integration in Eq.(35) has the lower and upper limits, ki ow and k upper , as the cutoffs of the 
wavenumber. The corresponding frequencies are vi ow ~ 2 x 10~ 18 Hz and v upper ~ 10 10 Hz. Detailed 
analyses of these limits are given in Ref.[26]. In past, in the absence of direct detection of RGWs, the 
constraint on RGWs through the energy density Qgw nas Deen commonly used; especially, a bound from 
the Big Bang nucleosythesis (BBN) [39] 

n GW h 2 < 8.9 x 10~ 6 (36) 

has been frequently employed in practice, where ho ~ 0.7 being the Hubble parameter [18]. 

In the following we give the resulting spectra h{k,Tn) and £l g {k) of RGWs, demonstrate their explicit 
dependence upon the model parameters 0, f3 s , 7, and the modifications by the QCD transition and e + e~ 
annihilation. Also we will compare it with the sensitivity of detections, such as LIGO and LISA. 

Figure 3 gives the spectrum h{v,TH) for three values of the inflationary index j3 = —1.8, —1.9, and 
-2.02, respectively, where the fixed r = 0.22, ft A = 0.75, (3 S = -0.3, n = 1.634, and v = 0.063 are 
taken. It is seen that h(u,Tn) is very sensitive to (3. A smaller (3 will generate lower amplitude of RGWs 
for all frequencies. 



Figure 4 shows the influence of reheating stage on RGWs. The spectrum is given for three different 
values of (3 S = 0.5, 0, and —0.3. It is seen that the reheating process will affect RGWs only in very high 
frequency range v > 10 7 Hz. This covers the frequency range of some very high frequency gravity wave 
detection systems, such as the Gaussian laser beam detector aiming at the frequency v ~ 10 10 Hz [40], 
or the circular waveguide detector aiming at v > 10 5 Hz [41]. We remark that the portion of predicted 
spectrum for the frequency range v > 10 10 Hz may be not reliable. This is because the energy scale of the 
conventional inflationary models are usually less than 10 16 Gev [26], which will gives a cutoff of frequency 
around v ~ 10 10 Hz. 

Presented in Fig. 5 is the modification of the spectrum h{v,Tu) by the QCD transition. There is a 
critical frequency vq ~ 10~ 9 Hz, below which the spectrum is the same for both models with or without 
the QCD transition. But, in the high frequency range v > u q Hz, the amplitude of h{v,Tn) is reduced 
by ~ 20% in comparison with the model without QCD transition. 

hk{T H )noQCD/h k (T H )QCD ^1-2 for V > Vq. (37) 

Note that the frequency range of this reduction on RGWs covers those of the major laser interferometers 
GW detectors, such as LIGO and VIRGO operating effectively in the frequency bands around ~ 10 2 Hz, 
and LISA around ~ 10 -3 Hz, respectively. So the modifications on RGWs due to the QCD transition 
are relevant to these major detections. In Appendix B we will give an interpretation to the origin of this 
critical frequency uq = 10~ 9 Hz and of the reduction in amplitude by ~ 20% for v > v q Hz. 

Fig. 6 gives the spectrum h(u,TH) modified by the e + e~ annihilation, which has a similar reduction 
effect. The lower frequency portion v < 10~ 12 Hz of the spectrum is not affected by the e + e~ annihilation. 
However, in the high frequency range v > 10~ 12 Hz, the amplitude of h{v,Tu) is reduced by ~ 10% by 
the e + e _ annihilation. Also the frequency range of this reduction covers that of LIGO, VIRGO, LISA, etc. 

Fig. 7 is an enlarged portion of the spectrum h{v,r\u) around v = 10~ 9 Hz, illustrating the details of 
reductions of the spectrum by both QCD transition and e + e~ annihilation. The combination of the QCD 
transition and the e + e~ annihilation reduce h(u,TH) ~ 30%. 

The influence of the dark energy on the spectrum h(y,Tn) is demonstrated in Fig. 8, where £1a = 0.0, 
0.7, and 0.75 are taken respectively. As explained in the last section, over the whole range of frequency 
(10 -19 ~ 10 10 ) Hz, the amplitude of spectrum is suppressed by the presence of S7a, but the slope remains 
the same. In particular, h(v,Tn) with Qa = 0.75 is reduced by ~ 45% in comparison with the model 
n A = 0. (See Eq.(29)) 

Figure 9 is a comparison of the LIGO detection with our calculated RGWs with the fixed tensor/scalar 
ratio r = 0.22 and the dark energy Q,\ = 0.75. The upper smooth curve is from the LIGO HI Upper limits 
(95% CL) from PowerFlux best-case [42], and the lower three fluctuating curves are the RGWs spectra of 
(3 = —1.8, —1.9, and —2.02, corresponding to those in Fig. 3, respectively. Here the vertical axis is the 
root mean square amplitude per root Hz, which equals to 

^. (38) 

The plot gives only the frequency range (30 ~ 300) Hz, on which the LIGO works efficiently. The 
rpdnrtinns Hup tr> the OTD transition and the annihilation havp hppn inmrnnratpd in thp rnrvp nf 



calculated h(u,Tn)/ 'yjv. Our result shows that there is a gap of about one order of magnitude even for 
the (3 = —1.8 inflationary model. As it currently stands, the possibility for LIGO to detect the RGWs 
predicted by the f3 = —1.8 inflationary model is not high, let alone other models with (3 < —1.8. Other 
two model parameters, i.e., the tensor/scalar ratio r and the dark energy SIa will substantially influence 
the height of h{u). The five-year WMAP data improve the upper limit on the tensor/scalar ratio r < 0.43 
(95% CL) for power-law inflationary models and r < 0.58 (95% CL) for models with a running index, 
and give the value of dark energy Q\ = 0.721 ± 0.015 [36]. So if we take the upper limit r = 0.43 for 
power-law models, the height of h(y) will increase by ~ 40% by Eq.(33), and, furthermore, if we take 
SIa = 0.721, h(u) will increase by another ~ 16% [3]. These together will allow an increase of h(v) by 
a total ~ 62%. Therefore, the current LIGO with greatly enhanced sensitivity [4] will definitely be able 
to put a constraint on the (3 = —1.8 inflationary model. However, it is seen from Fig. 9 that the curve 
for (3 ~ -2.02, supported by the scalar inflation models and the WMAP data [18] [36] [37], is about five 
orders below the LIGO sensitivity. So we may say that the RGWs generated by scalar inflation models is 
unlikely to be directly detected by LIGO at the moment. 

Figure 10 is a comparison of the LISA sensitivity curve with the spectra from Fig. 3 in the lower 
frequency range (10~ 7 , 10°) Hz that is also covered by the modifications of the QCD transition and the 
e + e~ annihilation. Assume that LISA has one year observation time corresponding to frequency bin 
Au = 3 x 10~ 18 Hz (i.e., one cycle/year) around each frequency. To make a comparison with the sensitity 
curve, we need to rescale the spectrum h(y) in Eq.(31) into the root mean square spectrum h(u,Au) in 
the band Au [31] [39], 



This r.m.s spectrum can be directly compared with the 1 year integration sensitivity curve that is down- 
loaded from LISA [43]. The plot shows that LISA by its present design will be able to detect the inflationary 
model of (3 = -1.8. If the ratio r > 0.22, it is still possible for LISA to detect the model of (3 = -1.9. 
However, LISA is unlikely to be able to directly detect the model of (3 = —2.02, as there is a gap of two 
orders. One may say that, in regards to detection of RGWs, LISA will give a stronger constraint on the 
RGWs spectrum than LIGO will do. 

Figure 11 shows the /3-dependence of the spectral energy density 0, g (u) defined in Eq.(34). Clearly, 
tt g is very sensitive to the inflationary parameter (3. A larger (3 gives a higher Q g . The Advanced LIGO [4] 
will be able to detect RGWs with U g h 2 > 10~ 9 at v ~ 100Hz, and it might impose stronger constraints 
on other inflationary models. On the other hand, the total energy density VLqw defined in Eq.(35) has 
also been used as a constraint on RGWs. Taking the parameters r = 0.22, J1a = 0.75, (3 S = —0.3, 
and n = 1.634, we find Qgw = 1-H x 10~ 2 for the inflationary model of (3 = —1.8, which is 3 orders 
higher than the BBN bound [39] in Eq.(36). So this model is disfavored, unless some other mechanism 
is introduced to reduce its £Igw- We also obtain £Igw = 1-97 x 10~ 8 for the model (3 = —1.9, and 
Q GW = 1.41 x 10~ 14 for the model (3 = -2.02; both models are safely below the BBN bound in Eq.(36). 

Fig. 12 shows £l g {v) in the (3 = —2.02 modified by the QCD transition and e + e~ annihilation, where 
the reduction on RGWs is more noticeable. In the range v > 10~ 9 Hz the QCD transition alone reduces 

0,.(iA hw ~ .30% anH thp rnmhinatinn nf OCR 1ran<;il-inn and thp p+p — annihilal-inn rpdnrpt; 0~(iA hv 




(39) 



~ 50%. 

The impact of dark energy on £l g {y) for the model [5 = —2.02 is plotted in Fig. 13. Three values of 
Q\ = 0.0, 0.70, 0.75 are taken. A larger gives a lower Q g . It is seen that the S7a causes a decrease 
of the amplitude of Q g by a factor of 1.6 over the whole range of frequencies. 
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Appendix A: Modelling of QCD transition and e + e~ annihilation 

The detail of the QCD transition is notoriously complicated and still under study. Here we will only 
consider the change of the effective degree of freedom, and give a simple working model for the scale 
factor a(r) around the QCD transition lasting a short period. Thus the radiation era can be tentatively 
divided into the three parts with the scale factor a(t) been listed in Eqs.(4), (5), and (6). After using the 
continuity conditions of a(r) and of a(r)' at the two given joining points t x and T2, one still needs to 
determine 3 parameters, r q , t x , n. Since the QCD transition temperature T oc 1/a during the radiation 
era, as soon as T is given, the initial time r q of QCD transition is determined, so is the parameter ( s . Still 
n and t x need to be fixed in the following. 

Around the QCD transition, the contributions from the matter p, m and the dark energy p^ can be 
neglected. Only the radiation is important, of which the energy density, the pressure, and the entropy 
density are given by 

2-7T 2 

Pr(T) = -^9*{T)T\ (40) 

Pr(T) = \ Pr {T), (41) 

s{T) = (p r +p r )/T = ^ S (T)T 3 , (42) 

respectively, where g* and g* s denote the effective number of relativistic species contributing to the energy 
density and entropy, respectively [44]. In an adiabatic universe, the entropy per unit comoving volume is 
conserved 

S(T) = s(T)a 3 (T) = constant. (43) 
From Eqs. (40), (42) and (43), one has 

p r oc g*g*s^a~ A . (44) 

During the QCD transition era, one has g* = g* s [29, 30]. Therefore, here we will not distinguish the 
difference between g* and g* s . Thus, Eq.(44) is reduced to 

p r oc (/* 1 ^ 3 a~ 4 , (45) 

which is plugged into the Friedmann equation (^-) 2 = yielding 

(a') 2 oc g~ 1/3 . (46) 



Table 1: The Values of n and £ 9 for different f- (T=190MeV) 





n 




1.3 


4.556 


1.25 


1.5 


1.634 


1.34 


1.8 


0.899 


1.47 


2.0 


0.714 


1.55 



Thus, with the decreasing of across the QCD transition the expansion rate a' in terms of the conformal 
time increases. As predicted by the Standard Model of particle physics, before the QCD transition, 
9*{Tq) = 51.25 at the time r q , and after the QCD transition, it becomes g^{r x ) = 17.25 at the time t x 
[45, 46]. Applying Eq.(46) to Eqs.(4) and (6) yields the ratio 

-9*(Tq)- 1 



a e 



1.2. (47) 



-9*(t x )- 

From the expression ( q = ^^y, one has 

n 

Cq 1+n =1.2. (48) 

If t x is determined (equivalent to giving Q), Eq.(48) will fix the parameter n. Let us determine ( q . The 
form of power-law in Eq.(5) implies g* oc t~ 6 ™ during the QCD transition era [30], thus Eq.(45) reduces 
to 

p oc T 2n a-\ (49) 
In terms of the cosmic time t defined by dt = adr, Eq.(5) implies 

l + n 

aocf 2 +». (50) 



Therefore, one gets 

c 9 = £ H = (r)^- ( 51 ) 

a{T q ) \t q y 

Once the ratio of t x /t q is given, then two parameters n and ( q can be fixed from Eqs.(48) and (51) 
immediately. Ref.[45] gives a value of the ratio t x /t q ~ 2; Ref.[46] estimated the typical duration of QCD 
transition to be the order of ~ 0.1t ? , i.e., t x /t q ~ 1.1. For concreteness in actual calculation, we take 

t x = (1.3 ~ 2.0)t q . (52) 

In Table 1 the values of n and £ q are listed for each given value of t x /t q . Here our treatment is different 
from that in Ref.[30], where the choice of n was put in by hand. 

In our model, there are three physical parameters: the QCD transition temperature T (i.e. the 
beginning time t q ), the duration of transition (i.e. the ratio t x /t q ), and the change of (i.e. the behavior 
of a(r) during the transition). As our calculation shows, £l g {y) is not very sensitive to the concrete values 
of T and of t x /t q . For T = 160 and 190 Mev, the difference in is only ~ 1%, and, similarly, for 

f.-/t- = 1 3 and 9 thp Hiffprpnrp in O..J1A k nnlv ~ 1% Thk <;mall Hiffprpnrp i<; qkn rpflprtpd hw thp tnt^l 



Table 2: ^gw for different QCD transition models 



with QCD transition 


£W(160M e y) 


n GW (190MeV) 


t-x/tq = 1-3 
tx/tq = 1-5 
W*g = 2.0 


1.40 x 10~ 14 
1.38 x 10~ 14 
1.35 x 10~ 14 


1.44 x 10~ 14 
1.41 x 10~ 14 
1.36 x 10~ 14 


without QCD trnsition 


1.54 x 10~ 14 


1.54 x 10~ 14 



energy density tt GW for different T and t x /t q listed in Table 2 for r = 0.22, fi A = 0.75, = -2.02, and 
Ps = —0.3. Therefore, as a conclusion, the most important element affecting the spectrum of RGWs is 
the change of g*, which are completely determined by the Standard Model of particle physics. 

The same treatment can be applied to the e + e~ annihilation. According to the Standard Model of 
particle physics, before the e + e" annihilation, g*(r y ) = g* s {T y ) = 10.75 at the time r y \ and after the 
e + e~ annihilation, g*(r z ) = 3.36 and g* s {r z ) = 3.91 at the time t z [29, 30]. Applying Eq.(44) to Eqs.(6) 
and (8) yields the ratio 

-1/2/ \ 2/3/ >. 
<h = 9* (Ty) g*s [Ty) _ 1 1 

a f 9 * 1/2 (r z ) gl' s \T z ) 

One can obtain 

Cj^ = LI- (54) 

Since T oc 1/a and the e + e~ annihilation occurred during T ~ (0.5 ~ 0.1) MeV [30], one obtains 
C y = ^ = 5 and v = 0.063 that are used in Eq.(7). 

Appendix B: Simple interpretation of modifications of RGWs by QCD transition and 

e + e~ annihilation 

Given the wave equation in Eq.(23), there are two limiting cases for wavelength. For the short 
wavelength k S> ^, the mode function h k has a decreasing amplitude, 

h k (r) oc l/a(r). (55) 

and, for the long wavelength k <C ^, the mode h k is simply a constant 

h k (T) ~ C k . (56) 

This property itself will be able to qualitatively account for why the QCD transition induces a decrease in 
the spectrum h(y^Tu) in Eq.(37). In Fig. 14 the scale factor a(r) is plotted for the QCD transition era 
around T ~ 190 Mev, and it is seen that 

a(T) QC D = a(T) noQ cD, t < 10~ 9 , (57) 
a{T)QCD -1.2, r>HT 9 . (58) 



a{T)noQCD 



This result, together with Eq.(55), explains the reduction of h(v,Tn) by ~ 20% by the QCD transition, 
as explicitly shown in Fig. 5. We need to explain why the reduction occurs for v > v q = 10~ 9 Hz. Since 
at any time r those modes with 

^ > ~TT ~ ~~ (59) 
a(T) t 

will decay as in Eq.(55), one sees that only those modes with 

k > k Q = — (60) 
T q 

will be reduced by the QCD transition. According to Eq.(18), in the present universe the physical frequency 
corresponding to kQ is 



® 2na(TH) 27T7r ( 



10" 9 Hz (61) 



for 7 = 1.044 and the Hubble constant Hq ~ 2.36 x 10~ 18 Hz. This critical frequency vq in Eq.(61) 
matches that demonstrated in Fig. 5 and Fig. 12, and is about 2 orders lower than 10~ 7 Hz, a value quoted 
in Ref. [29]. Thus, the features in the calculated spectrum of RGWs have been fully explained. We remark 
that a higher temperature T of the QCD transition would yield an earlier time r q and a higher critical 
frequency vq. As for the reduction by e + e~ annihilation, the interpretation is similar to the above. 
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Figure 1: The derivative o!{r) jumps up around the QCD transition modelled by Eq.(5). If there 
were no QCD transition , a!(r) would be a constant in this graph. 
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Figure 2: The scale factor a(t) in the accelerating, and non-accelerating models, respectively. Note 
that the horizontal axis is the cosmic time t. At the present time t = 1 the ratio of scale factors is 
aACDM/acDM — 1.8. 
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Figure 3: The spectrum }i{v,th) of GRWs is very sensitive to the index (3 of inflation. 
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Figure 4: The reheating affects the spectrum h{y^Tu) only in very high frequency range v > 10 7 Hz. 
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Figure 5: The QCD transition reduces h{v : T}j) by ~ 20% in the high frequency range v > 10 9 Hz, 
which covers the frequency band of LIGO and LISA. 




Figure 6: The e + e annihilation reduces h(u,TH) by ~ 10% in the frequency range v > 10 12 Hz, 
which covers the frequency band of LIGO and LISA. 
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Figure 7: The spectrum h(u, tjj) around v = 10 9 Hz for three cases: (1) with no QCD transition nor 
e + e~ annihilation, (2) with QCD transition only, (3) with both QCD transition and e + e~ annihilation. 
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Figure 8: The spectrum h(u,TH) depends upon the dark energy Q\ in the accelerating universe. A 
larger J7a yields a lower h(v,Tji)- 
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Figure 9: Comparison of the spectra with the LIGO HI Upper limits (95% CL) from PowerFlux 
best-case in the Ref.[42]. The vertical axis is the r.m.s amplitude per root Hz defined in Eq.(38). 
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Figure 10: Comparison of the spectra with the LISA sensitivity curve [43]. The vertical axis is the 
r.m.s spectrum defined in Eq.(39). The inflationary model of j3 = —1.8 can be tested by the LISA. 
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Figure 11: The spectral energy density Sl g {y) for various values of the parameter (3, 
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Figure 12: The spectral energy density Sl g (v) for three cases: (1) with no QCD transition nor e + e 
annihilation, (2) with QCD transition only, (3) with both QCD transition and e + e~ annihilation. 



12 



r=0.22 p^-0.3 
P = -2.02 n = 1.634 v = 0.063 




log i0 v [Hz] 



Figure 13: The spectral energy density £l g {v) for different A larger J1a yields a lower Q. g (y). 
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Figure 14: The scale factor a(r) is affected by the QCD transition. For r < 10~ 9 , cl(t)qcd = 
a{T~)noQCD- But for r > 10~ 9 , a(T H ) QC D / a{T H ) noQC D ^ 1-2. 



